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ABSTRACT 

Dynamical friction leads to an orbital decay of massive objects like young compact star clus- 
ters or Massive Black Holes in central regions of galaxies. The dynamical friction force can be 
well approximated by Chandrasekhar's standard formula, but recent investigations show, that 
corrections to the Coulomb logarithm are necessary. With a large set of N-body simulations 
we show that the improved formula for the Coulomb logarithm fits the orbital decay very well 
for circular and eccentric orbits. The local scale-length of the background density distribution 
serves as the maximum impact parameter for a wide range of power-law indices of -1. . . -5. 
For each type of code the numerical resolution must be compared to the effective minimum 
impact parameter in order to determine the Coulomb logarithm. We also quantify the cor- 
rection factors by using self-consistent velocity distribution functions instead of the standard 
Maxwellian often used. These factors enter directly the decay timescale and cover a range of 
0.5. . . 3 for typical orbits. The new Coulomb logarithm combined with self-consistent veloc- 
ity distribution functions in the Chandrasekhar formula provides a significant improvement of 
orbital decay times with correction up to one order of magnitude compared to the standard 
case. We suggest the general use of the improved formula in parameter studies as well as in 
special applications. 
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1 INTRODUCTION 

Dynamical Friction is a subject with two faces. Many researchers 
think it is s ufficiently well unders tood and studied since the classi- 
cal work o f lChandrasekhaiU 19421) and the many follow-ups. There 
is a duality between a collective gas-dynamic approach, study- 
ing responses in a continuum through which a test body moves, 
and a kinetic or particles approach where we consider a test par- 
ticle moving through a sea of light particles. As is k nown since 
long jBondi & Hoyl3ll944lRephaeli & Salpeterlll980h under cer- 
tain limits both approaches can yield similar results. In kinetic 
approaches (rooted in Chandrasekhar's work) usually an infinite 
homogeneous system of field stars is assumed, and the singular- 
ity at large impact parameters is cut off by the use of a so-called 
Coulomb logarithm In A (see more detail below). The underlying 
approximations are that interactions extend to impact parameters 
large compared to r g = GMbh/3<7 2 (Mbh the test mass particles, 
a the 1 -dimensional field particle velocity dispersion) and that the 
relative velocities involved are not too small. The error made by 
the first assumption is usually absorbed by fitting a certain nu- 
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merical value of the Coulomb logarithm to results of numerical 
simulatio ns; such a method has been very su ccessful in plasma 
physics (Rosenbluth, MacDonal d. Judd 1 19571) . star cluster dy- 
namics jSpitzerll 19871; iGiersz & Heggielll994h. galactic dynamics 
(iBinnev & TremainJl987[:IWeinberJl989l:lH ei-nquist & Weinberg! 
ll989l:|Prugniel & Combesll 19921 ; ICora. Muzzio. Vergnelll997h and 
in lSpinnato et alj J2003I) . where all citations should be seen as ex- 
emplary rather than exhaustive. 

There are two developments which have caused renewed inter- 
est in more accurate theoretical determinations of dynamical fric- 
tion. One is the realisation (both from theoretical structure for- 
mation models and HST observations of galaxy cores) that most 
cores of galaxies in the standard picture of hierarchical struc- 
ture formafiOT_are_embecMedj^ cuspy distribution s of dark mat- 
ter dNavarro. Frenk. Whitd[l997l : lLauer et alJI 19951) . It means that 
in many galaxies the density profile of the dark matter, which is 
the main background for dynamical friction of dwarf galaxies, star 
clusters and compact objects is nowhere constant as assumed in the 
standard Chandrasekhar theory. This problem has led to the sug- 
gestion of an empirical variation of the Coulomb logarithm with ra- 
dius, so_jtsJoac£c™t fOTdied^ of dynamical fric- 
tion jTremaindll976l : lHasfiimoto et alj|2003h . In ljust & Penarrubia 
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d2005l) a detailed theoretical investigation of the parameter depen- 
dence in Chandrasekhar's dynamical friction formula was given. 
An improved analytic approximation for the Coulomb logarithm 
was presented. 

Recently some doubt has been cast on the validity of 
Chandrasekhar's dynamical friction formula at all for certain 
densit y profiles, s uch a s harmonic potential, constant density 
cores. iRead et al] d2006h claim that while in other density pro- 
files it may work wel l (eve n in flattened systems, such as 
IPenarrubia. Just. Kroupal d2004l) have shown, if corrections due to 
the velocity distribution are taken into account), harmonic cores 
are supe r resonant structures where re sonances suppress dynamical 
friction. iBovl an-Kolchin et al. I d2008l) compare estimates of merger 
time scales of galaxies within dark matter halos obtained from 
Chandrasekhar's dynamical friction formula with results of high 
resolution iV-body simulations. They find that Chandrasekhar's 
friction formula does not work well, though they attribute it to fi- 
nite size and mass loss effects of the merging objects rather than 
resonant processes. Also llnouel |2009) supported the result that dy- 
namical friction ceases in harmonic cores, but doubt about resonant 
effects as the reason without giving further explanation. 

The idea that dynamical friction tends to stop in harmonic 
cores is a phenomenon not so uncommon in other fields of stel- 
lar dynamics. For example in simulations of super-massive black 
holes in galactic nuclei ([porband. Hemsendorf. & Merrittl [20031 ; 
iGualandris & Merrittl I200&). dynamical friction brings the black 
hole into the core, where it then starts a wandering process at- 
tributed to random encounters / Brownian motion. Also in star clus- 
ter dynamics dynamical friction of heavy mass stars sinking to the 
centre will stall if random encounters support enough pressure in 
the core. We think it is unclear whether it is a resonant or a ran- 
dom process which stops dynamical friction in cores, but this prob- 
lem deserves more detailed analytical and numerical attention. In 
numerical simulations it was shown dNakano & Makinolll999^lbh 
that in elliptical galaxies the settling of a SMBH to the centre can 
form a shallow cusp in the stellar distribution of the core. In these 
simulations the influence radius of the SMBH was not resolved and 
there are also very f ew observatio ns resolving the influence radius. 
The investigation of Merritt] d2006l) on the orbital decay of a mas- 
sive secondary BH points in the same direction creating a shallow 
inner cusp. 

In this paper, we assume Chandrasekhar's dynamical friction 
ansatz as a working hypothesis to be valid at least in principle, but 
after studying the velocity di s tribution dependence in earlier pa- 
pers dJust & Penarrubia 2005; Penarrubia, Just, Kroupa 20041) we 
focus here on the correction due to strong local density gradients, 
i.e. the opposite case to the disputed harmonic core situation. The 
detailed investigation of shallow cusps is postponed to a forthcom- 
ing paper. The generalised formula for dynamical friction is valid 
for extended objects like satellite galaxies or star cluster and for 
point-like objects like SMBHs. For extended objects the Coulomb 
logarithm is small and corrections to the relevant impact parame- 
ter regime are more significant than for SMBHs. Additionally the 
mass loss and the determination of the effective mass for dynamical 
friction must be taken into account dFuiii. Funa to. Makinol [20061 ; 
iFuiii et al]|2008l) . The present investigation is restricted to the or- 
bital evolution of SMBHs. 

Super-massive black holes, most likely to be present 
in merging galaxies fr om t he early universe o nwards 
dKormendv & Richstond 1 19951 ; iHaehnelt & Reesl 1 19931 ; 



iFerrarese et alj 2001), will sink to the centres of galac- 
tic merger remnants by dynamical friction and ultimately 



coalesce themselve s i Fukushige, Ebisuzaki. Makinol 1 19921 ; 

iMakino & Ebisuzakil Il99a ~ iBerentzen et alj 120091) . Numerical 
simulati ons to follow this process in a p a rticle-by-particle ap- 
proach dBerczik. Merritt. & Spurzen] 20051; iBerczik et al. 20061; 
Makino & Funatol |2004|; Hemsendorf. Sigurdsson. & Spurzeml 
120021 ; iMilosavlievic & Merrittll200lh are still too computationally 
expensive for realistic particle numbers, and so this situation 
requires another careful look at dynamical friction. Here not only 
the question of the influence of the ambient density gradient is 
important, but the test particle (super-massive black hole, single or 
binary) will typically violate also the other condition mentioned 
above, since it has small velocity (tendency towards equipa rtition) 
at leas t in the final phase of its approach to the centre. In Merritt 
d200ll) we find a comprehensive overview of dynamical friction for 
such objects in the limit of small velocities. An overview of how 
density gradients affect the sinking time-scales of massive black 
holes in galactic nuclei due to dynamical friction has apparently 
never been carried out. This is the aim of the present paper. 

Dynamical friction time-scales are also important for another 
aspect of massive black hole (binary) dynamics, namely the ec- 
centricity of the binary. Dynamical friction in homogeneous me- 
dia tends to circularise initially eccentric binary orbits, since it 
is most efficient at low velocity in apo-centre. Density gradients, 
if treated with stand ard (constant) Coulomb loga rithm should in- 
fluence this effect. Tsuchiva & Shimada (2000) provide a nice 
overview of how dynamical friction (including effects of den- 
sity gradients and anisotropic velocity distribution) affects orbital 
shap e (eccentricity), but they use the local epicyclic approxima- 
tion. lHashimoto et alj d2003h have shown that the circularisation 
problem can be solved by adopting a position dependent Coulomb 
logarithm. Recent models of three black holes in galactic nuclei 
show that dynamical friction together with three-body dynam- 
ics (e.g. Kozai effect) can induce e xtremely high eccentricities 
of th e innermost black hole bin ary dlwasawa. Funato. & M akino 
2006; Amaro-Seoane et al]|2010h . Also, the cosmological growth 
of massive central black holes from minor and major merging de- 
pends sensitively on dynamical friction of satellite galaxies and 
massive black holes in a background of stars and dark matte r 
dVolonteri. Haardt. Madaul2003l ; l\^lonteril2007l ; lDotti et al]201Ch . 

Direct numerical simulations resolving the full range of im- 
pact parameters and relative velocities for taking dynamical fric- 
tion correctly into account are still extremely tedious. Therefore it 
is very important to improve our theoretical ansatz for dy namical 
friction in non-standard cases. In I Just & Penarrubi3 d2005l) an im- 
proved analytic approximation for the Coulomb logarithm was pre- 
sented. In this article we concentrate on dynamical friction of com- 
pact objects in stellar cusps covering a wide range of parameters for 
circular and eccentric orbits. We present a comprehensive numer- 
ical analysis to test the general applicability of the Chandrasekhar 
formula with the new Coulomb logarithm. We take also into ac- 
count the self-consistent velocity distribution functions entering the 
friction force. 

In Sj2]the cusp models, the cumulative distribution functions 
and the Coulomb logarithm are discussed. In *j3]the theory of the or- 
bital decay is presented. In [J4]the iV-body and semi-analytic codes, 
which we are using to evolve our models, are described. In fj5] a 
comparison of the numerical results with semi-analytic calculations 
is given. In Sj6]a brief discussion of some applications is presented 
and finally Sj7]includes concluding remarks. 
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2 DYNAMICAL FRICTION FORCE 

Throughout the paper we use specific forces, i.e. accelerations. The 
dynamical friction force of a massive object with mass Mbh and 
orbital velocity Vbh in a sea of lighter particles is usually deter- 
mined by adopting locally a homogeneous background density p 
and an isotropic velocity distribution function. The result is a drag 
force anti-parall el to the motion, given by the standard formula of 
Chandrasekhar dBinnev & Tremainell 19871) 

Vdf = TJ2 X mA Wltn X= (1) 

In general the functions \ an d A depend on the velocity of the 
massive object and on the properties of the background system. 



2.1 Cusp models 

In order to quantify the position dependence of the Coulomb loga- 
rithm we investigate the orbital decay of a massive object in stellar 
cusps. We investigate two general scenarios. In the self-gravitating 
case the gravitational force is dominated by the power-law cusp. 
Bulges and the cores of galaxies and of star clusters can be de- 
scribed in that way as long as the massive object is outside the 
gravitational influence radius of a central black hole (or additional 
mass concentration). The circular speed is coupled to the cusp mass 
distribution by the Poisson equation. 

In the Kepler case the gravitational force is dominated by a 
centrally concentrated mass M c , which can be a central SMBH of 
the core of a self-gravitating stellar distribution. If the galactic nu- 
cleus already harbours a SMBH with mass M c at the centre, the in- 
ner part of the cusp inside the influence radius is dominated by the 
Kepler potential of the SMBH. In this case the stellar distribution 
can be described by a Bahcall-Wolf cusp. The orbital evolution of 
a second BH entering this region of gravitational influence is quite 
different to the self-gravitating case. This will happen, if a small 
galaxy merges with a larger one, both harbouring a central BH. In 
the minor merger process the bulges of the galaxies will relax first 
to the new bulge and the massive SMBH will settle to the centre. 
Finally the second BH decays to the centre by dynamical friction. 

Also in the outskirts of self-gravitating systems the density 
distribution may be approximated by a power law and the potential 
by a point-mass potential, if it is dominated by the mass concen- 
trated at the centre. We investigate two cases with steep power law 
distributions to test the maximum impact parameter dependence of 
the Coulomb logarithm. The Plummer sphere with an outer density 
slope of —5 and the Dehnen models with a slope of —4 are the most 
extreme cases. 



2.1.1 Power law profiles 

For analytic estimations we use idealised power law distributions 
for the background particles. For both cases (self-gravitating and 
Kepler potential) the density profile of the cusp and the cumulative 
mass profile can be approximated by the functions 
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It is comfortable to normalise all quantities to the initial values of 
the orbit. The initial position of the object is Ro = R(t — 0) 
leading to yo = 1, the enclosed mass is Mo = M(yo), and the 



circular velocity at Ro is V c , o = V c (Rq) = \JGMo/Rq. For the 
self-gravitating case we simply set M c = leading to Mo = Mt 
via Eq.[3] In the Kepler case we must distinguish between M (y) m 
M c = const, for the gravitational forces leading to the circular ve- 
locity with Vc 2 = GM c /Ro y~ x and the local density determined 
by Mt. The case rj = 5/4 corresponds to the Bahcall-Wolf (BW) 
cusp o f stationary equilibrium with constant radial mass and energy 
flow iBahcall & Wold 19761 : IQghtman & Shapiro! 19761) . 

We normalise all velocities to the local circular velocity V c 
instead of the velocity dispersion a by 



= v/V c and U = V hh /V c 



(4) 



The velocity dispersio n in self-gravitating c usps behaves quite dif- 
ferent for different rj iTre maine et alj |"l994). For r\ > 2 the kinetic 
pressure pa 2 converges to a finite value at the centre, which de- 
pends on the outer boundary conditions. The transition case with 
n = 2 is of special interest, because it corresponds to the p oc y~ x 
cusp as in the standard NFW cusp and in the Hernquist model. In 
this case the kinetic pressure pa 2 oc (C — lny), where the con- 
stant C is also determined by the outer boundary conditions. In a 
Kepler potential the isotropic distribution function degenerates for 
rj — 5/2, because it is completely dominated by particles with low 
binding energy. Dependent on the outer boundary conditions the 
distribution function can vary between a 5-function at the escape 
velocity (i.e. at E = 0) and a power law oc \E\~ with some cut- 
off at E rs 0. For a secondary black hole on a bound orbit in an 
idealised isotropic cusp the x- va l ue in Eq. [TJ becomes very small 
for rj > 2 and tends to zero for rj — 5/2 leading to an unreal- 
istically small dynamical friction force. Any change in the outer 
boundary conditions, small perturbations to the isotropy, or the ef- 
fect of higher order terms in the dynamical friction force become 
important. 

U can be converted to the standard X variable with the nor- 
malised circular velocity X c by 



X = X C U with x 2 c =^- 

Inserting Eq. |A9| into Eq.0 we find 



X 




rj < 2.5 Kepler 

rj < 2 self-grav. 

rj = 2 self-grav. 

2 < rj ^ 3 self-grav. 



(5) 



(6) 



which is independent of position y for the Kepler potential and the 
self-gravitating cusp with rj < 2. For a shallow self-gravitating 
cusp with rj > 2 the integration constant C depends on the outer 
boundary conditions of the realisation of the cusp. The details of 
the distribution functions are described in App.lAl 



2.1.2 Physical models 

In many simulations of stellar cusps it turned out that the setup of 
an initial cusp distribution in dynamical equilibrium with an un- 
physical outer cutoff is not stationary. The density profiles evolves 
deep into the inner cusp region. Therefore it is necessary to set up 
initial particle distributions and velocities with a well-defined outer 

cutoff of the cusp distribution. 

For the self-gravitating cusps we use Dehnen models toehnenl 

1 19931) with an outer power law slope of — 4 for the density . Thes e 
models are identical to the 77-models of ITremaine et al.l d 19941) . 
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Density and cumulative mass are given by 
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with Eq.[3]for the conversion for po - The Jaffe and Hernquist mod- 
els correspond to 77 = 1 and 77 = 2, respectively. Well inside the 
scale radius a the particles behave asymptotically like in the ide- 
alised power law distributions. 

For the Kepler potential case we investigate two different sce- 
narios. In the first case we investigate power law cusps in the vicin- 
ity of a central SMBH with mass M c , i.e. the BW cusp and a shal- 
lower Hernquist (He) cusp. There are no exact distribution func- 
tions known describing the Kepler potential part inside the influ- 
ence radius of the SMBH and the transition to a self-gravitating 

I ■ 1 J 1 

outer regime. iTremaine et al . ( 1994) generalised their 77 models by 

including the gravitational potential of a central SMBH and de- 
rived the power law distribution function (Eqs, IAl| and |A10^ well 
insid e the influence radius, whi ch is comparable to the scale radius 
a. InlMatsubavas hi et al.l d2007l) this approximation was adapted to 
a Plummer model instead of a 77-model. This model has two ad- 
vantages. Firstly the steeper slope in the outer part saves particles 
and computation time for simulations of BW cusps. Secondly the 
power law distribution function of the Plummer sphere is exact also 
in the inner part. Therefore realisations with smaller SMBH masses 
relative to the cusp mass are closer to equilibrium. The cumulative 
mass and density distribution are given by 
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where M t is the total mass of stars in the cusp. The approximate 
DF is 
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Here E is the binding energy and s = 5. We used these equations 
with 77 = 5/4 for the BW cusp. The constant fo represents the BW 
cusp, /1 the Plummer model and Eo is the transition threshold in 
energy. The He cusp is realised in a similar way. 

The outskirts of Dehnen models (DE) and the Plummer model 
(PL) can also be approximated by cusps in a Kepler potential using 
asymptotic expansions in y. From an identification of the density 
slopes for Dehnen and Plummer of -4 and -5, respectively, with 
770 — 3 in Eq. [3] we find 770 = —1 and 770 = —2 for the outskirts 
(here we use the index in order to distinguish it from the parame- 
ter 77 in the core). This leads with Eq. lAlOl to the correct distribution 
functions in Eq. lAll for the Dehnen models and the Plummer sphere 
in the limit of small energies. If we now identify Mt in Eq.[2]with 
the mass deficiency compared to the total mass Af c , then Eqs.|2]and 
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Figure 1. The plot shows the relative variation x(t/) / 'x B (U) as a function 
of 77 for for different values of U, where Xs corresponds to the standard 
Maxwellian distribution function (i.e. rj = 1). For self-gravitating cusps the 
full line is for the circular velocity [7 = 1, dotted and dot-dashed lines are 
for U = 0.7 and U = 1.4, respectively. The circles give the corresponding 
values for the Kepler potential cases with increasing size for increasing U 
(open symbols are for positive 77 and full symbols for negative 77). 



[3]also hold for these cases. We find for the Dehnen models 
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and similarly for the Plummer sphere 
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completely consistent with the power law cusp description. 



2.2 Cumulative distribution functions 

The cumulative distribution function x(U) of the normalised 1- 
dimensional distribution function F(u) measures the fraction of 
background particles with velocity smaller than U = Vbh/V c . We 
are using the self-consistent \ functions which are derived in ap- 
pendix [A] for the different models. They are significantly different 
to Xs(U) of the standard Maxwellian which is usually adopted for 
Chandrasekhar's formula. In most applications of the standard for- 
mula the local velocity dispersion is not known. Instead X c = 1 as 
in the singular isothermal sphere is adopted leading to the identifi- 
cation of X = U. In Fig.[TJthe correction factors \{U) /Xs(U) en- 
tering the friction force formula Eq.[Tjare shown. The different lines 
give the results for the self-gravitating cusps as a function of 77. We 
show the values for the circular velocity (7=1 (full line) and for 
U — 0.7, 1.4 (dotted and dot-dashed line), typical values for apo- 
and peri-centre velocities, respectively. In shallow self-gravitating 
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cusps (with rj 1) the efficiency of dynamical friction is reduced 
roughly by a factor of 77 due to the larger fraction of high velocity 
particles. In steep cusps dynamical friction is larger compared to 
the isothermal case, but with systematic deviations from the simple 
scaling for higher velocities U during peri-centre passage due to 
the finite escape velocity. The circles in Fig.[T]show x{U)/xs{U) 
for the BW and HE cusp (open symbols) and the outskirts of the 
Plummer (PL) and Dehnen (DE) spheres (full symbols) at the cor- 
responding values for rj — 5/4, 2, — 2, —1, respectively. For the 
HE case we used the numerically realised values (see also figure|4] 
For circular orbits the orbital decay time varies up to factor 
larger than two compared to the standard formula due to the self- 
consistent x functions. For the evolution of the orbital shape, the 
relative variation of the friction force between apo- and peri-centre 
is also important (see Sect. 15.1.31 ). 

2.3 Coulomb logarithm 

The main uncertainties in the magnitude and parameter dependence 
of the dynamical friction force is hidden in the Coulomb logarithm 
In A, which gives the effective range of relevant impact parameters 
and is up to now a weakly determined quantity. Since this formula 
is applied to wide ranges of parameters, it is very useful to have the 
explicit parameter dependen ce of In A instead of fittin g a constant 
value for each single orbit. In lJust & PenarrubiH |2005) the effect of 
the inhomogeneity of the background distribution on the dynamical 
friction force with Chandrasekhar's approach was discussed. The 
authors derived the approximation 
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(16) 



The Coulomb logarithm depends on the maximum and minimum 
impact parameter o max and o m in, resp., and for point-like objects on 
ago, the typ ical impact parame t er for a 90° -deflection in the 2-body 
encounters. Just & Penarrubli] J2005h found that the maximum im- 
pact parameter is given by the local scale-length D r determined by 
the density gradient, i.e. 

fe max = D x = JL- = JL ^2. (17) 
|Vp| 3-7? 

In an isothermal sphere b max is a factor of 2 smaller than the dis- 
tance R to the centre. In shallow cusps with 77 > 2 the local scale 
length D r exceeds the distance to the centre. In that c ase, the lo- 
cal sca le-length may be substituted by R (but see also iRead et all 
(2006) for additional suppression of dynamical friction in homoge- 
neous cores). 

For point-like objects like BHs the effective minimum impact 
parameter 090 is given by the value for a 90° -deflection using a 
typical velocity 7Jt yp for the 2-body encounters 
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GMv, 
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-R = 



3r e 



These values are a property of the numerical code and do not de- 
pend on the application. They are determined by numerical exper- 
iments (see also below). For extended objects like star clusters a 
good measure of the minimum impact parameter fe m i n is the half- 
mass radius TV 

The parameter dependence of b m&K and ago leads to a position 
dependence of In A, which affects the decay time Tdcc (see Eq.|29t 
and which also reduces the circularisation of the orbits resolving a 
longstanding discrepancy between numerical and analytical results. 
Using the distance to the centre a s maximum impact parameter was 
propo sed by different authors dTremainel 1 19761 : lHashimoto et al.l 
2003), but the effect on orbital evolution was never investigated 
in greater detail or for larger parameter sets. 

On circular orbits the local scale-length D T and the deflection 
parameter ago are position dependent. On eccentric orbits ago de- 
pends additionally on the velocity. Therefore In A varies systemat- 
ically during orbital decay and for eccentric orbits along each revo- 
lution. This has consequences on the decay time and on the evolu- 
tion of orbital shape. In eccentric orbits the dynamical friction force 
varies strongly between apo-and peri-galacticon mainly due to the 
density variation along the orbit. The variation due to higher peri- 
centre velocity and to the position dependent Coulomb logarithm 
weakens the differences. All these factors depend on the slope of 
the cusp density. That means that the effective Coulomb logarithm 
averaged over an orbit depends differently on the eccentricity for 
different values of rj. 

In the case of circular orbits with constant X c the position 
dependence of the Coulomb logarithm (eq.U6t can be parametrised 
by 

In A = ln(A /) . (20) 

Deep in shallow self-gravitating cusps with 77 > 2 the contri- 
bution from the circular velocity vanishes and A is also described 
by eq.|20] With Eq.[l8]the Coulomb logarithm is 



A 



extended or unresolved 



(3-ij) 6„ 

for extended objects and for point-like objects we find 



(21) 



/3 = 
P = V 
/3 = 4- 



T] 



A = 
A = 
A = 



(6-17) 



(3^,)(4-»)) M bh 



M, 



C" M, 



M h 



1L 



Kepler 
self-grav. , 
self-grav. , 



77 < 2 

2 < 77 < 3 



(22) 



We see that the motion of a point-like object in a Kepler potential is 
also described by a constant Coulomb logarithm, because the linear 
dependence of D r and ago cancel. The standard case corresponds 
to /3 = with Ro instead of D r and r g instead of ago in equation 
[16] leading to 



lnA s = 



In 



In 



Rg 



M 



unres. or ext. 



point-like. 



(23) 



, typ ^ T , bh .-r.. M(y)~ 2(1 + *, 

If the motion of a point-mass is numerically derived by a code, 
where ago is not resolved, the minimum impact parameter is deter- 
mined by the effective spatial resolution of the code. In our simula- 
tions we use the direct particle-particle code (PP) 0GRAPE with 
softening length e and the Particle-mesh code (PM) SUPERBOX 
with grid cell size d c . We use 



1.5 e (PPcode) 



d c /2 (PM code) 



( ) In a recent investigation ISpinnato et all d2003l) determined 

in a series of N-body calculations quantitatively the value of the 
Coulomb logarithm. They calculated the orbit of a massive point- 
like object moving through a background of stars with a power law 
cusp close to a singular isothermal sphere using different numer- 
ical codes (including SUPERBOX) and different parameters. They 
adopted a constant Coulomb logarithm and found that b max is sys- 
tematically about a factor of two smaller than the initial distance 
Ro to the centre. Additionally they restarted the orbital evolution at 
(19) smaller initial distances and find a smaller best-fitting In A. For ago 
they also found a value significantly smaller than the kinematically 
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defined gravitational infl uence radius r g of the moving BH. An in- 
spection of the results o f ISpinnato et~atl fe003l) shows, in spite of 
their slightly different interpretation, that they are fully consistent 
with equations [T6]-[22] /3 « 1 leading to a linear decrease of A 
for the resolved and the unresolved cases; & max = D T ~ Ro/2; 
femin = rfc/2 for the PM code SUPERBOX; a 90 = 0.75 r g . 



3 ORBITAL DECAY IN GALACTIC CENTRES 

We consider the orbital decay of a massive object in central stellar 
cusps in detail. We investigate the effect of the varying Coulomb 
logarithm In A and of a self-consistent distribution function f(E) 
instead of the standard Maxwellian on the decay rate. The dynam- 
ical friction force (Eqs.[]} depends on the background distribution 
via the local density of slow particles p(< Vbh) = PX» me density 
gradient, i.e. D r , and the circular velocity V c . 

Here we give the explicit solution for a massive object moving 
on a circular orbit in a power law density profile. The background 
distribution is a self-gravitating stellar cusp with the corresponding 
self-consistent distribution function f(E) or a stellar cusp in a Ke- 
pler potential. We use xo = x(U = 1) for the circular speed. The 
initial decay timescale for angular momentum loss is given by 



to 



Vco 
Vdf.o 



Ml 



27rr;xo mA M bh M, 



To 



(24) 
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Ro Mo' 
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_pc M Q _ 




M M 



Myr. 



The first expression is in units of the orbital time To = 27r7?o/V c o 
and the second expression is in physical units. In the case of a self- 
gravitating cusp the enclosed cusp mass Mo equals M t and to is 

1/2 

proportional to M t . In the Kepler case Mo equals M c and ro is 
proportional to M c . 

For the radial evolution y(t) of circular orbits with varying 
Coulomb logarithm (0 ^ 0) we find the implicit solution (Ap- 
pendix |Bj 



t = 



2 k 



r A1 — \Ei(]nz )-Ei(\nz(y))\ k,P^0 
zo 

K A k/B 

z y = A >' 
I 3 — 2rj Kepler 



(25) 



3 + 7] self-grav. 
and for the special cases the explicit solutions 

k^0,/3 = 



»(*) 



1 - -E- 



A (cxp(-i/ Tdf )-l)//3 K = 0;/3 ^ 



Kepler 



(26) 



k = P = Kepler 
The angular momentum evolution can be easily calculated by 



L(t) = a/ GM(y)R = La x | ^ 



(1+,))/2 self-grav. 
1 Kepler 



(27) 



with I/o = \J GMoRo- We introduced 



3-2r; 



self-grav. 
Kepler 



Tdf = To X < 



K = o, p Kepler 



(28) 



k = P = 



Kepler 



and | Tdf | measures the decay timesca le. The first line of eq .l26lre- 
produces the special solution given in lSpinnato et alj d2003l) . 

For shallow cusps with r\ > 3/2 in a Kepler potential k be- 
comes negative leading to a negative Tdf. In that case In zq is also 
negative and there is formally a stalling of the orbital decay for 
P > when A approaches unity. In case of P = equation [26] 
yields an infinite decay time to the centre. 

Also for positive k the total decay time with varying In A 
is not well-defined, because the approximations in eq. 1 161 breaks 
down for In A < 1. In order to get an analytical estimate of the 
effective decay time Td cc , we choose the time needed to decrease 
In A from the initial value InAo to In A = 0.3725/3/k, where 
Ei(\n z) = 0. For In zo 3> 1 it can be estimated from eq. 1281 with 
the help of the asymptotic expans ion Ei(x) ~ e x (l/x + l/x 2 ) 
(8.216) dGradshtevn. Rvzhik 1983) 



Tde 



1 + 



K In Ao 



T d f for 



InAo > 1 . 



(29) 



The correction factor quantifies the effect due to the position de- 
pendence of the Coulomb logarithm, if P ^ 0. For a negative k 
(as realised in the Hernquist cusp HE) the last term in eq. [^dom- 
inates and diverges as A approaches unity. Therefore we define a 
decay timescale for a fixed minimum value j/dcc (e.g. three times 
the stalling radius) by 

lnz 



Trie 



ZO 



■ In I In z\ Tdt for I In z\ = 



lnA«l. (30) 



For k — and (5/0 the orbital decay of the BH would also 
stall at In A = 0. Only for k > and p = there is a finite 
time Tdf to reach the centre. For completeness we mention that for 
a negative p the total decay time to the centre would be finite due 
to the enhanced friction force by an increasing Coulomb logarithm. 
But for all realistic cases we find P ^ 0. 

The standard case of dynamical friction corresponds to the 
isothermal sphere and a constant Coulomb logarithm, i.e. initial en- 
closed mass Mo at radius Ro with f] = 1 and P = in Eq.[28] We 
use the decay time of the standard case 



TdfO = 



M Ro 



XslnA s Mbh Vco 



M<1 

Mt 



self-grav. 
Kepler 



(31) 



as given also in lBinnev & Tremaind dl987h (7-26) as normalisation. 
For more details of the orbital evolution see App.lBl 

3.1 Self gravitating cusps 

In Fig.[2]the orbital decay y(t) of circular orbits is presented. The 
differences in the orbital evolution are caused by the combination of 
using self-consistent density profiles and distribution functions and 
by the position dependence of the Coulomb logarithm. The stan- 
dard case with In A = const. (P — 0, r\ — 1) is given by the black 
broken-dotted line. For point-like objects y(t) is given in the top 
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Figure 2. Radial evolution of circular orbits in cusps. We have used 
M hh = 10~ 2 M and r h = 6.7 X 10~ 3 i?,o leading to the same In A = 
5.0 for the standard case with /3 = 0, 7] = 1. For the four Kepler cases 
Bahcall-Wolf cusp (BW), Hernquist cusp (HE), Dehnen (DE) and Plummer 
(PL) outskirts with r\ = 1.25, 2.0, —1.0, —2.0 we have chosen the initial 
cusp masses Mt/M c = 0.5, 0.5, —0.225, —0.06, respectively. The top 
panel is for point-like objects and the lower panel for extended objects. 
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Figure 3. Relative variation of total decay times (eq. |29t normalised to the 
standard case of constant Coulomb logarithm. We used the same parameters 
as in Fig. |2] Open symbols are for positive rj and full symbols for negative 
V- 



panel. Orbits for different power law indices r\ with varying In A 
are plotted. The full red line shows in an isothermal core the de- 
lay due to the position dependence of In A and the slightly smaller 
initial value of In Ao = 4.6 instead of 5.0. In the bottom panel the 
evolution for the same values of rj are shown for extended bodies. 
The parameters are chosen to give the same Coulomb logarithm 
In Ao = 5.0 in the standard case. 

In Fig. [3] the variation of the effective decay time as a func- 
tion of rj is shown for point-like and extended bodies. Some care 
should be taken to use these numbers, because the innermost radius 
reached at time Tdcc depends on rj. But the general parameter de- 
pendence of Tdcc gives some insight in the physics of the orbital 
decay in cusps. The effective decay time of circular orbits in self- 
gravitating cusps is affected by the following aspects: 

• Density profile With decreasing rj the mass Mo is more and 
more concentrated to the centre leading to a smaller density in the 
outer regions. This results in a smaller dynamical friction force and 
prolonged T dcc . 

• Distribution function With decreasing 77 the fraction of slow 
particles \o increases considerably from 0.2 to 0.7, reducing the 
effect of the smaller density in the outer parts. 

• Coulomb logarithm The position dependence of In A leads 
to a moderate delay in orbital decay in the later phase. The effect is 
strongest for large values of rj. 

• Extended bodies For extended bodies like star clusters In A is 
generally smaller compared to point-like bodies, because the min- 
imum impact parameter is larger. Compared to the standard case, 
the prolongation factor is only weakly dependent on rj. 

In the case of eccentric orbits there is an additional effect of 
the variation of In A, because along these orbits the relative strength 
of the friction force at apo- and peri-centre is changed. 



3.2 Kepler potential 

In case of a Kepler potential the cusp mass distribution is decou- 
pled from the potential and the enclosed cusp mass M t is an ad- 
ditional free parameter. We discuss the explicit solutions for four 
cases. Well inside the influence radius of a central SMBH the stel- 
lar distribution can be described by a cusp in a Kepler potential. 
We present the orbital decay in the Bahcall-Wolf cusp and the shal- 
low cusp of a Hernquist profile. In the outskirts of self-gravitating 
systems the density distribution may be approximated by a power 
law and the potential by a point-mass potential, if the density pro- 
file is steep enough. We investigate two cases with steep power law 
distributions to test the maximum impact parameter dependence of 
the Coulomb logarithm. The Plummer sphere with an outer density 
slope of —5 and the Dehnen models with a slope of —4 are the ideal 
cases. 



3.2.1 Bahcall-Wolf cusp 

Here we look to the orbital decay inside the gravitational influence 
radius, where the enclosed mass of the stellar component M t is 
smaller than the mass M c of the central SMBH. We neglect in this 
region the contribution of the stellar component to the mean gravi- 
tational field. 

The general equations are already given in the previous sec- 
tions, but we evaluate the terms explicitly for the Bahcall-Wolf cusp 
with 



rj = 5/4 p oc y 



-7/4 



(32) 



8 A. Just et al. 



leading to 
1 

K = — 

4 



8 



(33) 



For circular orbits the minimum impact parameter for a point-like 
object becomes 

11 Mbh 
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leading to the initial Coulomb logarithm 
P — point-like 
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The distribution function F(u) and the corresponding x{U) are 
shown in Figs. IA1I - 1 A3 1 The decay time-scale Tdf from Eq. [28l 
reads 
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(36) 



For point-like objects the total decay time Tdcc equals Tdf and for 
extended bodies the corresponding equation is (using Eq.|29t 
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lnA 



(InAo) 
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unres. or ext. 



The orbital decay of circular orbits is given by Eqs.l26land[ 
point-like and extended bodies, respectively, 

4 
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Tdf 



(37) 
5] for 

(38) 



which is very different to the standard case. 

In Fig. [2] the orbit evolution in a Kepler potential is plotted 
for M t = 0.5 M c . It shows the strong slow-down in the inner part 
leading to long total decay times Tdcc (see Fig. [3}. 



3.2.2 Hernquist cusp 

For the Hernquist cusp the corresponding equations are 
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The decay timescales are 
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and the orbital decay of circular orbits is for point-like objects 
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(44) 



In Fig.[2]we have chosen the same j/h as for the BW case. We find 
an initially faster decay than in the BW cusp and a stronger slow- 
down in the late phase leading to a comparable Td cc adopting a final 
Coulomb logarithm In A = 1.1 (see figure 0. The comparison of 
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Figure 4. Initial and final cumulative distribution functions for the Hern- 
quist cusp at different radii of runs El and E2. For comparison the analytic 
limit for r — ¥ is included. 



the analytic and the numerically realised cumulative distribution 
function x(U) is shown in Fig. [4] The figure demonstrates the in- 
fluence of the outer boundary conditions deep into the cusp. 



3.2.3 The outskirts of a Plummer sphere 

In the outskirts of a Plummer sphere with total mass M c the corre- 
sponding equations are 
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We have used eq.[l5]for Mt. The orbital decay of circular orbits is 
for point-like objects 
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and the distribution function F(u) and the corresponding x{U) are 
shown in Figs. lAll - IA3l In Figs.|2]and[3]we have chosen y a = 0.2 
for the Plummer radius in units of Ro. 



3.2.4 The outskirts of Dehnen models 

In the outskirts of Dehnen models with total mass M c and inner 
cusp slope r\ the corresponding equations are 

-4 



5 

K= 2 



pocy 
2 = 5 
2 

5 Mbh 



X. 



U9 ° = 714^' 



(52) 
(53) 

(54) 



Dynamical friction of massive objects 9 



An = 



7_ Mc 

20 A/ bh 



1 go 

4 *>™in 



/3 = point-like 
P = 1 unres. or ext. 



Tdf = 



Td 



2.36 
2/a In A 
lnA + 0.4 



'Ro~ 


3/2 


' M c " 


1/2 


'M bh " 


_pc_ 




_M Q _ 




_M _ 



"Tdf 



unres. or ext. 



(55) 



Myr, (56) 



(57) 



(InAo) 

We have used eq. [14] for M% with r\ — 1.5. The orbital decay of 
circular orbits is for point-like objects 
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and the distribution function F(u) and the corresponding x(U) are 
shown in Figs. lAll - IA3l In Figs.[2]and[3]we have chosen y a = 0.15 
for the scale radius in units of Ro . 



4 CODE DESCRIPTION 

We are using two different numerical codes to evolve our initial 
models. In order to study the orbital decay of the massive BH in 
self-gravitating cusps and in the outskirts of the Plummer sphere 
and Dehnen models, we used the PM code SUPERBOX. In order to 
study the decay of a secondary massive BH in a Bahcall-Wolf cusp 
around a SMBH, we used the PP code <f)GRAPE. For comparison 
of the numerical results with the analytic predictions for eccentric 
runs, the semi-analytic code INTGC has been developed. A brief 
description of these codes is given below. 

4.1 SUPERBOX 

The PM code SUPERBOX dFellhauer et al.| [2000) is a highly effi- 
cient code with fixed time step for galaxy dynamics, where more 
than 10 million particles per galaxy in co-moving nested grids for 
high spatial resolution at the galaxy centres can be simulated. The 
code is intrinsically collision-less, which is necessary for long-term 
simulations of galaxies. Another advantage of a PM code is the 
large particle number which can be simulated in a reasonable com- 
puting time (up to a few days or a week per simulation). Three grid 
levels with different resolutions are used which resolve the core 
of each component/galaxy, the major part of the component/galaxy 
and the whole simulation area. The spatial resolution is determined 
by the number of grid cells per dimension iV c = 2 m and the size 
of the grids. The SMBH is included as a moving particle with own 
sub-grids for high spatial resolution in its vicinity. The grid sizes 
are chosen such that the full orbit of the secondary BH falls into 
the middle grid. 

All SUPERBOX runs were performed using the Astronomis- 
ches Rechen-Institut (ARI) fast computer facilities with no special 
hardware. 

4.2 ^GRAPE 

This programme is a direct iV-body code, which calculates the pair- 
wise forces between all particles. In order to combine high numer- 
ical accuracy and fast calculation of the gravitational interactions 
a very sophisticated numerical scheme and a special hardware is 
used. This PP code uses fourth order Hermite Integrator with in- 
dividual block time steps. The acceleration and its time deriva- 
tive are calculated with parallel use of GRAPE6A cards. In order 



to minimise communication among different nodes MPI paralleli- 
sation strategy is employed. For the simulations in a dense cusp 
around a massive central SMBH, we are using a specially devel- 
oped 0GRAPE code. We include the central SMBH as an external 
potential in order to avoid the relatively large random motion of 
a live SMBH due to the small particle number. Secondly the time 
step criterion is modified. We add a reduction factor for the BH 
time step, which compensates the effect of the relatively small ac- 
celerations compared to the field particles. The cod e and special 
GRAPE hardware are described in lHarfstetal.ll2007h . For our cal- 
culations we are using the parallel version of the programme on 32 
node cluster Tifarfl built at the Astronomisches Rechen-Institut 
in Heidelberg. 

Part of the calculations were done on the special 85 node Tesla 
CI 060 GPU cluster installed on the National Astronomical Obser- 
vatories of China, Chinese Academy of Sciences (NAOC/CAS). 
For these runs we used the modified version of our <^GRAPE 
code including the SAPPORO library t o work on the GPU cards 
dGaburov. Harfst. Portegies Zw art 2009). 



4.3 Semi-analytic code - INTGC 

The program INTGC is an integrator for orbits in an analytic back- 
ground potential of a galactic centre including the Chandrasekhar 
formula for dynamical friction. Different analytic mod els with their 
X fun ctions and a variable Coulomb logarithm jjust & Penarrubia 
1 20051) are implemented. An 8th-order c omposition scheme is used 
for the orbi t integ ration jYoshidal 1 1990h : for the coefficients see 
iMcLachlanl ( 119950 ). Since the symplectic composition schemes 
are by construction suited for Hamiltonian systems, the dissipa- 
tive friction force requires special considerati on. It is implemented 
in IN TGC with an implicit midpoint method ( Mikkol a & Aarsethl 
l2002h . Four iterations turned out to guarantee an excellent accu- 
racy of the scheme. Gravitational potential and density are given 
analytically. 



5 NUMERICAL MODELS AND RESULTS 

The contributions to the dynamical friction force covers a large 
range of parameters for the 2-body encounters, which must be fully 
covered by the numerical simulations in order to reach a quantita- 
tive measure of the Coulomb logarithm. The numerical represen- 
tation is mainly restricted by the resolution of small impact pa- 
rameters determined by the spatial resolution, the time resolution 
and the number statistics. For PP codes the number statistics in the 
main limitation. Therefore we can use the <^GRAPE code for the 
BW cusp only with a very high local density in the inner cusp. We 
reach a few encounters per decay timescale to with impact param- 
eters below twice the minimum value. For the SUPERBOX runs the 
spatial resolution is limited, which we can take into account by a 
correct choice of the minimum impact parameter. But even in that 
case it turns out that due to the fixed time step the time resolution is 
still a bottleneck, which limits the total time of some simulations. 
In all models we use the angular momentum to measure the orbital 
decay (eq.|27|l. 



http://www.ari.uni-heidelberg.de/grace/ 
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Table 1. Parameters of the runs for the BW cusp. 



Run 


AT/10 3 


e/10" 4 




Ro 


V c ,o 


Vo/V c>0 


ago/10- 4 


InAo 


AO 


64 


1.0 












orm 


Al 


64 


0.1 


.005 


0.2 


2.25 


1.0 


5.79 


-§-.28 


A2 


64 


0.1 


.005 


0.1 


3.17 


1.0 


2.89 


.1.29 


A3 


64 


0.1 


.01 


0.2 


2.25 


1.0 


11.6 


1.59 


Bl 


64 


0.1 


.005 


0.2 


2.25 


0.7 


7.01 


&09 


B2 


64 


0.1 


.005 


0.1 


3.17 


0.5 


4.07 


S.94 


CI 


64 


1.0 


.005 


0.2 


2.25 


1.0 


5.79 


|.25 


C2 


64 


5.0 


.005 


0.2 


2.25 


1.0 


5.79 


■4.79 


C3 


64 


10.0 


.005 


0.2 


2.25 


1.0 


5.79 


4.26 


El 


64 


0.1 


.005 


0.7 


1.20 


1.0 


17.5 


5 33 


E2 


128 


0.1 


.002 


0.2 


2.27 


1.0 


2.0 


6.91 



5.1 Bahcall-Wolf cusp 

We used the extended 77-model of iMatsubavashi et alj J2007h to 
generate the initial particle distributions in phase space. In all 
our runs we are using N = 64, 000 particles. The particles po- 
sitions are generated so that their spatial distribution satisfy Eq. 
(|9) and the velocities were assigned to these particles according 
to Eq. (TO). In all our simulations we are using the normalisation 
G — M c = a = 1 leading to y a = Rq 1 in N-body units. The mass 
of the central black hole is M c and the total cusp mass M t = 0.1. 
For the setup we used a radius range of 10 -4 — 20 in all our runs. 
Inside one unit length our cusp follows density profile —7/4 and 
then turns over to a Plummer density profile with slope —5 for far 
out distances. Table[T]shows the list of parameters for the series of 
runs. 



E 

Z5 

O 



0.1 



0.01 



0.001 



Note. For all runs the gravitational constant G, the SMBH mass M c and the 
scale radius a are normalised to unity. N is the total number of particles 
in the cusp, e the softening length of the particles, Afth tne mass of the 
secondary black hole with initial distance Ro in units of a, the circular 
velocity V c ,o at Ro in units of yj GM c /a and the initial velocity Vo in 
units of the circular velocity, the initial value of ago from Eq.[l8]is in units 
of a. 
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Figure 5. Top: Figure shows the evolution of Lagrange radii of 0.1, 0.5, 1, 
3, 5, 10, 30, 50, 80, 90% enclosed mass (from bottom to top) for the BW 
cusp. The Lagrange radii do not show any systematic evolution with time. 
Bottom: Cumulative mass profile at various time steps. The cumulative 
mass profile is practically indistinguishable from the theoretical one except 
the inner few dozen particles, where deviations due to the inner cutoff and 
noise are expected. 



5.1.1 Cusp stability analysis 

Since we are using an approximate DF, we first performed a run 
(run AO in the table) without a secondary BH to test whether or 
not the cusp is stationary around the central SMBH. We run this 
model up till 50 time units. Figure [5] shows the evolution of La- 
grange radii and also the cumulative mass profile at various time 
steps. We can see that the cusp is very stable. So the problem how 
to get a stationary cusp in a Kepler potential was solved by this 
kind of a compromise distribution between Bahcall-Wolf cusp in 
the inner (high binding energy) regime and a Plummer distribution 
in the outskirts. 

5. 1.2 Circular Rims 

We performed two series of circular runs (see table [TJ. The runs 
A1-A3 with different BH masses and initial radii are all resolved, 
i.e. the softening parameter e = 10~ 5 is much smaller than the 
initial value of ago. Since ago decreases linearly in y (Eq.|34|>, the 
minimum impact parameter is fully resolved to very small radii 
and we can use Eq.[38]for analytic estimates. In figures |6]and|7] we 
show the distance and the angular momentum evolution, because 
already small deviations from circularity smear out the appearance 
of the orbits due to the very short orbital time compared to the 



decay time. The top panel of figure |6] shows the distance evolution 
for runs Al, A2, and A3. Run A2 corresponds to a restart of Al 
after time T = 155 but using the initial particle distribution of 
the cusp. We see that there is no significant long-term evolution of 
the cusp, which influences the orbital decay, until the end of run 
Al. The dashed (green) lines show the analytic predictions of the 
orbital evolution from Eq. [38] There is an excellent agreement in 
the first phase with a small delay in the later phase of run Al, which 
occurs much earlier in A2. The reason for the reduced friction is 
investigated in run A3 further. In run A3 we increased the mass 
of the BH and put it back at Ro — 0.2 such that in case A3 the 
radius, where the enclosed cusp mass Mt equals Mbh (horizontal 
lines in figures|6]and|7](, is twice that of run A2. In run A3 the delay 
starts also very early. In the bottom panel of figure [6] and in figure 
[7] the same evolution is shown much clearer in angular momentum 
L using Eq. 1271 for the analytic predictions. The horizontal lines 
show the distance, where Mbh = Mt in all figures. This radius 
coincides with the distance, where the BH mass exceeds the mass 
in a shell centred at the orbit. An inspection of the cumulative mass 
profiles shows that the back-reaction of the scattering events to the 
cusp distribution becomes significant (see figure[8}. The cumulative 
mass profile becomes shallower, which means that the local density 
is reduced leading to a smaller dynamical friction force. 
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Figure 6. Top panel: Comparison of the orbit evolution of the </>GRAPE 
data and the analytic estimates (Eq. |38t for the circular orbits Al, A2, A3. 
The run A2 is shifted by To = 155 in order to continue the theoretical line 
of Al. The horizontal lines shows the radii, where the enclosed cusp mass 
equals M bh (for Al, A2 at T > 150 and for A3 at T < 100). Bottom 
panel: Comparison of the angular momentum evolution of the </)GRAPE 
data and the analytic estimates (Eq. l27t for the circular orbit Al. The hor- 
izontal line shows L c at the radius, where the enclosed cusp mass equals 
M hh . 



For a circular orbit in a Bahcall-Wolf cusp the new friction 
formula is very close to the standard formula: the Coulomb loga- 
rithm is also constant, the value deviates only by A In A=-0. 1 from 
In A s , the \ value is 0.430 instead of Xs=0.428. This leads to an 
indistinguishably faster decay when applying the standard formula. 
The picture changes slightly for the eccentric orbits (see below). 

In all simulations the eccentricity of the orbits vary slightly. 
The increasing eccentricities in the later phases of the runs may 
correlate to the decreasing mean density, i.e. may be connected to 
the feedback of the BH on the cusp. The eccentricity evolution will 
be investigated in more detail in future work. 

In a second series of runs C1-C3, we study the impact of the 
minimum impact parameter by increasing the softening parameter 
until it is much larger than ago. In figure [9] we can see that for the 
largest softening length e (greater than (290) the decay is slower as 
expected from the smaller Coulomb logarithm (eq.|35l>. We can use 
the variation to measure the effective resolution of the code. The 
best simultaneous fit of all curves lead to 6 m i n = 1.5e in Eq.|16|for 
the PP code 0GRAPE. The random variations of the orbital decay 
on time-scales of 10. . . 100 time units are expected from a rough 
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Figure 7. Same as in the bottom panel of figure[6]for runs A2 and A3 
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Figure 8. Cumulative mass profiles for the runs A1-A3 demonstrating the 
reduced local density in runs A2 and A3 at the position of the BH. 

estimate of the close encounter rate and the corresponding velocity 
changes. 

5.1.3 Eccentric runs 

In eccentric runs the additional effect of the velocity dependence of 
x{U) and also 090 (eq.U8b along the orbit occurs. We performed 
two runs Bl and B2 with different initial radii and velocities at 
apo-centre corresponding to eccentricities of e = 0.5, 0.75, re- 
spectively. The effective minimum impact parameter ago is well 
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Figure 9. Evolution of a circular orbit in Bahcall-Wolf cusp for with differ- 
ent softening lengths e (runs Al, C1-C3). Thin lines show the semi-analytic 
calculations with INTGC with eq. |16| for the Coulomb logarithm with the 
best fit minimum impact parameter £> m i n = 1.5e. 

resolved along the orbits. In the top panel of figure [10] the orbit of 
B 1 is shown in two short time intervals of length AT = 1 in or- 
der to resolve individual revolutions. Since the orbital phase is very 
sensitive to the exact enclosed mass and apo- peri-centre positions, 
the cumulative phase shift after hundreds of orbits is expected. In 
the bottom panel the evolution in L is shown for runs B 1 and B2 
for the full evolution time. The comparison with the semi-analytic 
results show very good agreement for run B 1 . The delay in the or- 
bital evolution as in runs Al -A3 does not occur. In run B2 there 
is a delay in the later phase. In contrast to run B 1 the peri-centre 
passages of B2 suffer from low number statistics. The peri-centre 
distance decays from 0.015 at T = to 0.01 at T = 150, where 
only 200 cusp particles are enclosed inside the orbit. 

In the lower panel of figure[l0]we show also the orbital evolu- 
tion using the standard parameters. In the standard case with con- 
stant In A s (equation|23t and Xs from a Gaussian velocity distribu- 
tion are applied. The decay is slightly faster. In order to separate the 
effect of the new Coulomb logarithm and the correct distribution 
function we show also the orbital decay substituting only x by Xs 
(dotted blue line) and only In A by In A s (dashed-dotted cyan line), 
respectively. The effect of the new Coulomb logarithm is larger, but 
still not significant in the Bahcall-Wolf cusp. 

5.2 Hernquist cusp 

We performed two circular runs El and E2 in the shallow Hernquist 
cusp to test the limit of our friction formula using the <^GRAPE 
code. Run El is performed on the GRAPE cluster at ARI and E2 
with larger N on the GPU cluster of NAOC/CAS. The orbital decay 
is shown in figures[TT]and[T2] The velocity distribution function de- 
viates significantly from the analytic limiting case of an idealised 
cusp (see figure |4). The \ function is stable over the simulation 
time and depends only weakly on the distance to the central SMBH. 
Therefore we used for the semi-analytic simulations with INTGC 
constant mean values for the circular orbits. The local scale length 
of run El entering the Coulomb logarithm is smaller than the dis- 
tance to the centre, which is the limiting value, because Ro = 0.7 
is close to the scale radius. The orbital decays are very well repro- 
duced in R and in L. A comparison with the standard values In A s 
and Xs show a significant deviation mainly due to the different \ 
function. 
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Figure 10. Top Panel: Orbit evolution of run B 1 in two time intervals (0, 1 ) 
and and (100,101) compared to the semi-analytic predictions from INTGC. 
Bottom panel: Angular momentum decay for runs B 1 (top) and B2 (bottom) 
for the full simulated time compared to the semi-analytic predictions. Label 
INTGC stands for the new In A and \, 'standard' for In A s and Xs, the other 
two for substituting only one parameter by the standard value. 

5.3 Outskirts of the Plummer sphere 

The background distribution is a self-gravitating Plummer sphere 
and the outskirts are described by eq. [15] We place the orbits of 
the circular run PI and eccentric run P2 far outside the Plummer 
radius in order to be very close to the power law density with 770 — 
3 = —5 and a Kepler potential according to eqs. [15] Since the 
number density is very small, we cannot use a PP code. Instead we 
use SUPERB OX with ten million particles. The parameters of all 
SUPERBOX runs are listed in table[2] 

Due to the very long orbital decay time the density distribu- 
tion evolves slowly. In steep cusps the representation of the local 
density in the semi-analytic calculations is the most critical param- 
eter. This affects the decay time scale via the local density and also 
the enclosed mass (eq.l24b. In the Kepler limit the enclosed mass is 
given by the constant central mass M c . In lowest order the density 
evolution can be modelled by a linear expansion of the Plummer 
radius. The analytic solution of circular orbits is still given by eqs. 
[25] and [26] with the substitution t — > s(t) with v = 2 (see eq IBlOb . 
By fitting the density profiles at different times we find t a — 11.3 
Gyr. 

The numerical results of the circular run PI are compared in 
figure [T3] to the analytic solution of eq.[25]in the Kepler limit us- 
ing an expanding Plummer profile. For the calculations with INTGC 
we correct additionally for the decreasing enclosed mass. Since in 
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Figure 11. Orbital decay of a BH in the Hernquist cusp of the circular orbit 
El in r(t) (top panel) and L(t)/Lo (bottom panel). The 0GRAPE simula- 
tion is compared to the semi-analytic results with INTGC. Same notation as 
in figure [10] 



these simulations ago is not resolved by the grid cell size d c , we 
determined the best value for the minimum impact parameter to be 
6min = do/2. The same factor 1/2 is used for all other SUPER- 
BOX runs. In the lower panel of figure \\3\ the orbits with standard 
Coulomb logarithm In A s and standard \s as in figurefTOlare shown 
for comparison. In figure [14] the orbital evolution is shown for the 
eccentric run P2. 

In both cases we find a satisfying agreement of the numerical 
results and our analytic predictions. With the standard formula the 
decay is significantly delayed (figures[T3land[T4t. Here the correc- 
tions due to the correct \ function and the new Coulomb logarithm 
have a different sign and cancel each other partly. In the circular run 
both corrections are equally important. In the eccentric run the cor- 
rection by the \ function dominates, but the additional correction 
due to the new In A is also significant. 



5.4 Outskirts of Dehnen models 

The background distribution is a self-gravitating Dehnen model 
with 77 = 1.5 in the inner cusp. We place the orbits of the circular 
run Dl and the eccentric run D2 (tablefj} far outside the scale radius 
in order to be very close to the power law density with 770 — 3 = — 4 
and a Kepler potential according to eqs.l 141 

The numerical results are compared in fi gures [73] and 1 1 6 1 with 
the semi-analytic calculations using INTGC. In the numerical simu- 
lations the inner cusp becomes shallower and shrinks considerably 
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Figure 12. Same as in figure fTTI but for run E2 with higher resolution and 
closer in. 

leading to a decreasing enclosed mass and increasing density in 
the outer parts. We correct for that in INTGC. In these simulations 
ago is not resolved as in the Plummer case. The numerical results 
are in good agreement with the analytic predictions using the same 
minimum impact parameter & m j n = d c /2. The orbital resonances 
seen in figure Q3] are connected to a motion of the density centre 
of the cusp, which is used as the origin of the coordinate system. 
Therefore the angular momentum of the orbit oscillates due to the 
accelerated zero point. 

5.5 Self gravitating cusps 

Here we study the orbital decay of a BH in self-gravitating cusps 
of Dehnen models using SUPERBOX without a central SMBH. We 
performed two runs (one circular and one eccentric) for a Hern- 
quist (7/ = 2) and a Dehnen-1.5 (77 = 1.5) model each. For the 
semi-analytic calculations we use again the same minimum impact 
parameter b m i n = d c /2 in Eq. 1 161 

5.5.1 Dehnen-1.5 model 

We studied the decay of a circular and an eccentric orbit in the in- 
ner cusp of a Dehnen model with r\ — 1.5 (D3 and D4 in table[2](. 
A comparison of d c /2 and ago shows that the close encounters are 
marginally resolved. For this reason and because the transition of 
the inner and outer power law regimes is very wide, the analytic 
approximation of a self-gravitating cusp shows systematic devia- 
tions. With INTGC we find nevertheless a very good match to the 
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Table 2. Parameters of the SUPERBOX runs. 



Run Mtot M a r cut N m dt d c M bh Ro V Ct o Vb/Vc.o ago InAo 

10 9 Mq 10 9 Mg kpc kpc 10 7 Myr pc 10 6 M kpc km/s pc 



PI 


1.0 


0.970 


0.1 


10 


I 7 


0.3 


16.1 


1.0 


0.7 


77.2 


1.0 


0.538 


2.87 


P2 


1.0 


0.970 


0.1 


10 


I 7 


0.3 


16.1 


1.0 


0.7 


77.2 


0.7 


0.869 


2.87 


DI 


1.0 


0.826 


0.1 


1 


I 6 


0.1 


16.7 


1.0 


0.4 


94.2 


1.0 


0.186 


2.48 


D2 


1.0 


0.826 


0.1 


1 


I 6 


0.1 


16.7 


1.0 


0.4 


94.2 


0.7 


0.207 


2.48 


D3 


1.0 


0.128 


1.0 


10 


I 8 


0.1 


2.78 


1.0 


0.3 


42.8 


1.0 


1.24 


4.75 


D4 


1.0 


0.128 


1.0 


10 


I 8 


0.3 


2.78 


1.0 


0.3 


42.8 


0.7 


1.70 


4.44 


HI 


1.0 


0.197 


1.0 


10 


I 7 


1.0 


16.1 


1.0 


0.7 


34.8 


1.0 


6.0 


3.68 


H2 


1.0 


0.197 


1.0 


10 


I 7 


0.3 


16.1 


1.0 


0.7 


34.8 


0.56 


9.7 


3.46 


H3 


1.0 


0.052 


1.0 


10 


I 8 


0.1 


4.76 


1.0 


0.3 


30.3 


1.0 


2.13 


3.64 


H4 


1.0 


0.052 


1.0 


10 


I 8 


0.1 


4.76 


1.0 


0.3 


30.3 


0.9 


2.39 


3.53 



Note. All quantities are given in physical units. Mtot is the total mass inside the cutoff radius r cu t , Mo is the enclosed mass at the initial distance Ro , a is 
the scale radius, N is the particle number, 2 m is the number of cells per dimension in each grid, and d c is the cell size of the middle grid. The initial value of 
090 from Eq. ll8l is calculated for the exact model. 
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Figure 13. Orbital decay of a BH in the outskirts of a self-gravitating Plum- 
mer sphere for the SUPERBOX runs of the circular orbit PI in r(t) (top 
panel) and in L(t)/Lo (bottom panel). The SUPERBOX run is compared to 
the semi-analytic results with INTGC. Same notation as in figure fTol The top 
panel shows additionally the analytic approximation from equation l25l 



SUPERBOX results of the circular run D3 (see Fig. 1 1 5t and for the 
eccentric run D4 (see Fig. 1 1 6b . In the upper panel of Fig.ll5lwe 
show also the circular run with a larger time-step in SUPERBOX, 
where for the close encounters with impact parameter comparable 
to the cell length d c the motion of the perturber are not resolved 
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Figure 14. Same as in figure[l3]but for the eccentric orbit P2. The top pan- 
els show r(i) for two time intervals with higher resolution and the bottom 
panel shows L z (t) for the full calculation. Same notation as in figure [Tol 
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Figure 15. Same as in figure [13] but for the orbit Dl in the outskirts of a 
Dehnen model. Same notation as in figure [TOl 



in time. This leads to a larger effective minimum impact parameter 
and a slower orbital decay. 



5.5.2 Hernquist model 

The Hernquist model with shallower cusp (rj = 2) is even more 
complicated, because is not constant and the x function de- 
pends on the outer boundary conditions. For the two runs HI (cir- 
cular) and H2 (eccentric) the parameter range is similar to the 
Dehnen- 1.5 case with marginally resolved ago. The orbits of H3 
(circular) and H4 (eccentric) are further in and the grid resolution 
is higher leading to a fully resolved ago. All orbits are reproduced 
reasonably well by INTGC using the co rrect x function and taking 
the correct velocity dispersion in ago (Tremain e et al.|[l~994l) into 
account (see figures [T9l-I22t. 



5.6 Velocity distribution functions 

The local velocity distribution functions are crucial for the dynam- 
ical friction force. Therefore we check here the numerical reali- 
sation of the distribution functions. The best way to measure the 
distribution function is to determine the x function in a spherical 
shell at the distance of the BH for different times. In figurel23lthe 
X functions are shown for a few examples. The top panel shows the 
Kepler cases and the bottom panel the self-gravitating cusps. 

In the Kepler case we present the circular runs in the Bahcall- 
Wolf cusp (BW, Al, A3), the Plummer case PI, and the Dehnen 
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Figure 16. Same as in figure [T5l but for the eccentric orbit D2. Same nota- 
tion as in figure [Tol 



case Dl. In the Bahcall-Wolf cusp the final distribution function of 
Al is very close to the theoretical line. The x functions of the ec- 
centric runs Bl, B2 are very similar and therefore not plotted here. 
In the runs A2 and A3, where the density profile flattens slightly, 
a bump in the x function at velocities below the circular velocity 
can be observed. It is more pronounced in run A3. The x value at 
the circular speed is not influenced dramatically showing that the 
orbital delay is caused entirely by the reduced local density. The 
distribution functions in the outskirts of the Plummer and Dehnen 
cases are very stable and well represented by the theoretical expec- 
tations. 

The lower panel of figure[23"lshows the x functions of the self- 
gravitating cusps for the circular runs HI and D3. Here we added 
the case of the eccentric run D4 to demonstrate the stability of the 
velocity distribution functions independent of the shape of the orbit. 

Overall the velocity distribution functions are very robust and 
well reproduced in the numerical simulations. 
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Figure 17. Orbital decay in distance (top panel) and angular momentum 
(bottom panel) for the circular orbit D3 in the self-gravitating inner cusp of 
the Dehnen-1.5 model. Same notation as in figure fTol In the top panel there 
is additionally the orbit with insufficient time resolution shown. 



6 APPLICATIONS 
6.1 The Galactic Centre 

The central region of the Galactic Centre can be modelled by a cusp 
with T\ = 1.2 and enclosed mass M = M(R ) = 1 ■ 1O 9 M at 
R = 200 pc, and a cen tral BH with mass M c = 2.6 • 1O 6 M 
dGenzel & Towneslll987l) . The gravitational influence radius is at 
R — 1.4 pc, where we may assume a slightly shallower cusp with 
rj — 1.25. If the central BH would have entered the cusp on a cir- 
cular orbit at some later time, the decay time from Ro = 200 pc 
to the centre would be HOMyr (Eq. [29] with /3 = rj). An interme- 
diate mass BH with A/bh = 1 ■ 10 4 Mq can reach the centre from 
Ro ~ 60 pc in 2.5 Gyr. For the final decay inside the influence 
radius of R = 1.4 pc the decay time is 15 Myr. 

In the central cusp of the Galaxy there are young star clus- 
ters like the Arches and the Quintuplet cluster with a projected dis- 
tance from the Galactic centre of about 30 pc and an age of a few 
Myr. The stellar mass is ~ 1 ■ iq 4 M with a half- mass radius 
of r h « 0.2 pc dCotera et al.lll996l ; iFiger e"tai]|l999l) . Assuming 
that the initial distance is Ro = 30pc we find with Eq. [22] for the 
initial Coulomb logarithm In Ao = 4.4 leading to the decay time 
i"dcc = 720 Myr. If we assume that the star clusters are still em- 
bedded in their parent molecular cloud with mass 1 • 1O 6 M and 
initial half-mass radius r^o ~ 3 pc, the decay time becomes con- 
siderably smaller. Adopting cluster mass loss linear in time we get 
Tdcc ~ 30 Myr, which is still large compared to the actual age of 
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Figure 18. Same as in figure[T7]but for the eccentric orbit D4. 



the clusters. This is in contrast to Ger hard d200lh . who investigated 
the infall of massive clusters, also embedded in giant molecular 
clouds, to the Galactic centre. He found much shorter time-scales, 
because he used unrealistically high values for the Coulomb loga- 
rithm In A w 10 - 20. 

6.2 Minor merger 

One important application is the orbital decay of the SMBHs af- 
ter the merger process of two galaxies. Lets assume a 10:1 merger 
with primary BH of mass M c = 1 • 1O 8 M and the secondary 
BH with Mbh = 1 ■ 1O 7 M . After violent relaxation of the stellar 
components and settling of the primary BH to the centre we adopt 
a shallow new central cusp with r\ — 1.75 at radii large compared 
to the gravitational influence radius of the central SMBH. With an 
enclosed mass M = M(R ) = 1 ■ 1O 1O M at Ro = lkpc the 
circular velocity is V c o = 207km/s corresponding to a velocity 
dispersion of ao = 293km/s (eq. [6f. For the secondary BH we 
find with Eq. [29] a decay time of 830 Myr. For the inner 500 pc 
with enclosed mass of M (500pc) = 3 ■ 1O 9 M the BH needs 
190 Myr. After reaching the gravitational influence radius of the 
central SMBH at 7? = 72 pc, the final decay takes 15 Myr (with 
Ed. 1361). If we compare these decay times with an isothermal model 
with the same enclosed mass at Ro = 1 kpc and adopting a con- 
stant Coulomb logarithm of In A = 6.9, we find the correspond- 
ing times 792; 200; 4.1 Myr, respectively. The total decay time is 
comparable, but in the final phase the decay is significantly slower 
in the Kepler potential compared to the standard isothermal approx- 
imation. 
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Figure 19. Same as in Fig. ll7l for the circular orbit HI in the self-gravitating 
inner cusp of the Hernquist model. 



7 CONCLUSIONS 

Classical dynamical friction as derived by Chandrasekhar is still of- 
ten used for important astrophysical applications, such as infalling 
small galaxies in dark matter halos or globular cluster systems 
around galaxies. It is also important to understand the co-evolution 
of super-massive black holes in galactic nuclei (after mergers) with 
their host galaxies. These recent new interests in dynamical friction 
problems have caused several new studies of dynamical friction in 
higher order than originally by Chandrasekhar (e.g. regarding the 
velocity distribution) or in a different physical framework (collec- 
tive modes in stellar or gaseous systems rather than particle-particle 
interactions). 

Chandrasekhar's dynamical friction formula has been proved 
in many ways, and agrees with the study of collective modes, but 
with notable exceptions, one which is being discussed in the litera- 
ture now is the possible absence of dynamical friction in harmonic 
cores, relevant for cosmological structure formation, avoiding too 
strong infall of small galaxies. Still in many papers today the clas- 
sical Chandrasekhar formula is applied by using a local isothermal 
approximation for the kinematics of the background system (i.e. 
the x function) and by fitting the Coulomb logarithm In A for each 
orbit. 

In this article we tested quantitatively the effect of using the 
self-consistent velocity distributi on function in y and a ge neral ana- 
lytic formula for In A derived bv ljust & Penarrubial d2005l) . We per- 
formed high-resolution numerical simulations, with particle-mesh 
and direct A-body codes, of the orbital evolution of a massive 
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Figure 20. Same as in Fig. ll9l for the eccentric orbit H2. 



black hole in a variety of stellar distributions. We investigated cir- 
cular and eccentric orbits in self-gravitating cusps (Dehnen models) 
and in cusps in a central Kepler potential (so-called Bahcall-Wolf 
cusps) and in the Kepler limit of steep power law density profiles in 
the outskirts of stellar systems. The background distributions cover 
a large range of power law indices between — 1 ... — 5 of the den- 
sity profile. 

The application of the self-consistent \ functions lead to cor- 
rection factors in the orbital decay times in the range between 
0.5 ... 3 (Fig. [TJ. The main new feature in the improved general 
form of the position and velocity dependent Coulomb logarithm 
(Eq. |16t is the local scale length D T of the density profile as max- 
imum impact parameter. In most applications the effect of the new 
In A is a significant delay in the orbital decay. A detailed compari- 
son with orbital decay using the standard values In A s and \s shows 
that the corrections are very different in the different cases leading 
to a general improvement of the orbit approximation. In a few cases 
like circular orbits in the Bahcall-Wolf cusp with a resolved mini- 
mum impact parameter the standard formula can still be used. But 
already for eccentric orbits a measurable difference occurs. We like 
to point out the generality of the new formula such that a fit of in- 
dividual orbits is no longer necessary. We find a general agreement 
in the orbital decay of the numerical simulations with the analytic 
predictions at the 10% level. This holds for circular as well as ec- 
centric orbits and in all background distributions in self- gravitating 
and in Kepler potentials. 

Another more technical finding concerns the best choice of the 
minimum impact parameter 6 m i n measuring the numerical resolu- 
tion. It should be reminded, that from the structure of In A, there is 




implicitly one common scaling factor for o max , o m i n , 190 free. So 
the normalisation of one of t hese quantities mus t be fix ed in order 
to determine the other two. In Just & Penarrubia (2005) we argued 
for ago and 6 max . Here we fixed 6 max = D T and proved that ago 
is the correct effective minimum impact parameter in the numer- 
ically resolved cases. On that basis we determined the numerical 
resolution for the different codes in terms of the softening length e 
or grid cell size d c . We find 6 m j n = 1.5e for the PP code consistent 
with the results of other authors and b m i n = d c /2 for the PM code 
SUPERBOX. 

Our result is very important for any conclusions about the 
rates of massive black hole binary mergers and black hole ejections 
in the course of galax y formation in the hierarch i cal structure for- 



mation picture (cf. e .g. Volonteri. Haardt. Madaul j2003h : IVolonteril 
d2007h : lDotti et alj 120 id) ). It helps also to interpret results of 
recent numerical investigations of the pr oblem of multiple black 
holes in de nse nuclei of galaxies dMilosavlievic & Merrittl 
2001; iHemsendorf. Sigurdsson7& Spurzem |2002| ; 



Milosavljevic & Merritt 2003; Berczik, Merritt, & Spurzem 2005; 
Berczi k et al.l 120061 ; iBerentzen et alj 120091 ; lAmaro-Seoane et al.l 
20101) . We did not investigate in detail the feedback of the de- 
caying BH on the background distribution. A general flattening 
of the central cusp is expected for mass ratios of order unity 
dNakano & Makinolll999al| bl; lMerritti r2006). The investigation of 
dynamical friction in shallow cusps is much more complicated, 
because the maximum impact parameter is not well defined and 
the outer boundary conditions influence the velocity distribution 
function deep into the cusp. 

The new formula can be used for extensive parameter studies 



of orbital decay, where a high numerical resolution is not possible 
to resolve the dynamical friction force. Even if we did not test the 
formula for extended objects like satellite galaxies or star clusters 
in this article, it should work equally well in these cases. Since In A 
is generally smaller in these applications, the correction due to the 
maximum impact parameter would be more significant. Additional 
effects li ke mass loss of the satellite galaxies must be taken into 
account dFuiii. Funato. Makindl2006l ; iFuiii et alj|2008l) . The new 
Coulomb logarithm works also fo r non-isotropic background distr i- 
butions, as was shown already in Pen arrubia. Just, Kroupal 12004). 
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Figure 23. Top panel: Initial and final x functions for the circular runs in 
the Bahcall-Wolf cusp (BW, Al, A3), the Dehnen case Dl, and the Plummer 
case PI. Bottom panel: Same for the self-gravitating cusps with Dehnen-1.5 
(D3, D4) and Hernquist profiles (HI). 
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APPENDIX A: DISTRIBUTION FUNCTIONS 

The distribution function of a power law cusp (eqs.[2]and[3} is given 
by 



f(E) = 



K\E\ P 

K\E\ P 

Kexp{-E/a 2 ) 



Tj < 2.5 

o s; 77 < 3 

77 = 1 



Kepler 

self-grav. (Al) 
self-grav. 



with normalisation constant K and energy E = $ + v 2 /2, where 
the zero point of the potential is at the centre for self-gravitating 
cusps with 77 > 1 and at infinity for the Kepler case and self- 
grav. cusps with 77 < 1. The dependence of the power law index 
p on 77 is different for the Kepler and the self-gravitating case (see 
below). We consider the two cases where $ is given by the self- 
gravitating potential of the cusp or the Kepler case with <J> dom- 
inated by the central mass M c . The natural normalisation of the 
velocity is ^/|2$|, which corresponds to the escape velocity for 
vanishing potential at infinity. For practical use it is more comfort- 
able to normalise the velocities to the circular velocity (it from Eq. 

We introduce the normalised 1 -dimensional distribution func- 
tion F(u) and the cumulative function x(U) by 



47T7J 2 /(£ , )d7J = 

X(U) = 



pF(u)du 

u 

F(u)du . 



which can be written in the form 

K'u 2 (l + Zu 2 Y 
K'u 2 /1 ' t>*V 



F{u) = 



l + £u 
K'u 2 exp (- 



77 < 2.5 
< 77 < 3 
77= 1 



(A2) 
(A3) 



Kepler 

self-grav. (A4) 
self-grav. 



with 



V 2 

s 2$ 



Since K' and £ are constant, F(u) and x{U) are independent of 
position y. 

The velocity dispersion can be obtained by integrating the 
Jeans equation in volving the second momen t of the velocity dis- 
tribution function (Bin nev & Tremainell 19871) 
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Figure Al. Normalised 1-dim distribution function F(u) for different 
cases as a function of normalised velocity u = v/V c (see Eq. IA4t . The 
sequence with r) = 0.5, 1.0, 1.25, 1.95 shows a decreasing maximum. The 
full line is the Maxwellian (rj = 1). The Kepler potential cases are labelled 
by PL for Plummer, DE for Dehnen outskirts, BW for the Bahcall-Wolf 
cusp, HE for Hernquist cusp. 



where the integration constant C" depends on the inner and outer 
boundary conditions. 



Al Self-gravitating cusps 

The constants in Eqs. lAll and lA4l are given by 



and 



V = 



3 + 77 
2(1-77) 



(A6) 
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77 = 1 
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Here we used the Beta function B(x,y) — T(x)T(y)/T(x + y) 
(see iGradshtevn. Rvzhikl dl980h (8.38)). In Fig. \M\ the distri- 
bution functions F(u) are shown for different values of 77 = 
0.5, 1.0, 1.25, 1.95 (with decreasing maximum). The thick line is 
the standard Maxwellian. For 77 < 1 the energy range is finite with 
u 2 < whereas for 77 ^ 1 the potential is infinitely deep 

allowing for all velocities. 

For the dynamical friction force we need the cumulative func- 
tion x{U), which gives the fraction of stars with velocity smaller 
than U. In Fig. IA2I x(U) is shown for the same set of 77-values 
as in Fig. I A 1 1 In order to see more clearly the difference to the 
Maxwellian distribution function the ratios x{U)/Xs(U) are plot- 
ted in Fig. [A3] For 77 < 1 the \~ values are larger and for shallower 
cusps with 77 > 1 the values are smaller. In Fig. lA4l we show x(U) 



Dynamical friction of massive objects 21 




• /// ..Kepler: PLti=-2.0 

DEti=-1.0 

BWti=1.25 

HE n=2.0 

self-grav. r|=0.5 

7/ J n=i .0 

n=i.25 

11=1.95 



2.5 



3.5 



0.8 



0.6 



0.4 



0.2 




0.5 



1.5 



In I 



Figure A2. x(^) f° r me same values of rj as in Kg. En 



Figure A5. Same as in Fig.|A4]but for fixed X = 0.7, f .0, 1.4. 
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Figure A3. Same functions as in Fig. IA2I but here normalised to the 
Maxwellian: x[U)/Xs{U) (increasing rj from top to bottom). 




Figure A4. The functions x{U) as a function of r) are shown for fixed 
U = 0.7, 1.0, 1.4 (Dashed blue, full red, dotted pink line, respectively). 
The circles give the corresponding values for the Kepler potential cases 
with increasing size for increasing U (rj = —1 for DE, ij = — 2 for PL, 
r\ = 1.25 for BW, r) = 2.0 for HE). Open circles are for positive r) and full 
circles for negative r\. 



as a function of rj for fixed U.U = 1 corresponds to the circular ve- 
locity, and U = 0.7, 1.4 are typical values for peri- and apo-centre 
velocities with moderate eccentricities. Fig. |A5| gives the same but 
for fixed X = 0.7, 1.0, 1.4 showing the systematic difference due 
to the different ^-dependence of X. 

The velocity dispersion is needed for X c (eq. [6j used in ago 
(eq.|18ll. Substituting M(y) and p(y) from eq.|2]and eq.[3]yields 



c (y) = 



GM n-i _ V£ 
2R (2- V )y 2(2-7)) 



I C"'y 3 - V 



0<r/<2 

r, = 2 (A9) 

2 < r) < 3 



where the constants C and C'" depend on the outer boundary con- 
ditions. 



A2 Kepler potential 

In the Kepler potential of a central SMBH (without the mean field 
contribution of the stellar cusp) we find for the constants 

1 3 

2 ; P =2 



-v- 



(A10) 



and 
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K' 



4nV2B(3/2, 1 + p) 

2 | C |3/2 
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(All) 



(A12) 



B(3/2,l + p) 

In a pioneering work IPeeblesI jl972l) analysed the structure of a 
stellar cusp in a Kepler potential of a central SMBH, which is sta- 
tionary for times large compared to the relaxation time. Unfortu- 
nately he derived incorrec t values for p and thus r\. The correct 
derivation can be found in Lightman & Shapiro (1976) using scal- 
ing arguments and in lBahcall & W olf ( 1976) using a Fokker-Planck 
analysis. The resulting so-called Bahcall-Wolf cusp (BW) is given 
by V — 1/4 leading to r\ — 5/4 and the well-known density profile 
p oc y~ 7 ^ 4 . In Figs. lAll - |A5| the corresponding functions or values 
are shown (dashed lines or full circles, respectively). In a Hernquist 
cusp (HE) with rj — 2 we have a shallow density profile p oc y" 1 
and findp = — 1/2 leading to a diverging 1-D distribution function 
at the escape velocity. The outskirts of Dehnen (DE) and Plummer 
(PL) distributions with densities p oc y~ 4 and p oc y~ 5 correspond 
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to 77 — — 1 and r\ = —2 with p — 5/2 and p = 7/2, respec- 
tively. The corresponding distribution functions and x values are 
also plotted in Fies.lAll-lA5l 



APPENDIX B: ORBITAL DECAY 

Here we derive explicitly the orbital decay of a massive object on a 
circular orbit in a cuspy density distribution with a position depen- 
dent Coulomb logarithm. All values are normalised to the values at 
the initial distance Rq. The distance to the centre is y — R/Ro- 
The orbital decay can be computed by identifying the angular mo- 
mentum loss due to the dynamical friction force (eq.Q~|l 



V : 



y 

TO 



M ln(A / 
M(y) InAo 



(Bl) 



dy 

where we have used the definition of the decay timescale to (eq. 
124b . the position dependence of the density (eq. [3}, replaced the 
square of the circular velocity by GM(y) / Roy with initial value 
GMq/Ro, and the parametrisation of the Coulomb logarithm (eq. 
120b . The enclosed mass M(y) on the right hand side and the circu- 
lar velocity V c (y) on the left hand side are different for the Kepler 
case and the self-gravitating case. 

Bl Kepler potential 

In the Kepler potential with constant enclosed mass M(y) = 
M = M c and 

Vc = V^y- 1 



we write eq. lBll in the form 
2 In (A02A 



r In Ao 
We define for 77 ^ 3/2 



A? =A V 



'/ 



and find after a little mathematical manipulation 

dz Ik zq 
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For the Kepler potential the results can be easily generalised 
to expanding or contracting cusps. If the density varies proportional 
to a power of time, i.e. p(y) = p(y,t = 0)(1 + t/t & ) v , then the 
differential equation |B3I with time dependent to can be converted 
back to the original form with initial po in to by 



dt 



ds ~ V ds 



2 In (Ao/) 
T In A 
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The implicit solution of the differential equation |B 3 1 is still given 
by Eq. lB6l with the substitution t — > s. For a Plummer sphere with 
linear increasing Plummer radius j/ a = J/ a o(l + t/t & ) we find for 
the outskirts v — 2 (eq.|15b. 

B2 Self-gravitating cusps 

In a self-gravitating cusp we have j3 = 77, 1, for a point-like 
object, an extended object, and a constant Coulomb logarithm, re- 
spectively (Eq.l20b. For the self-gravitating case we insert M(y) = 
Moy v and 



Vc = V: o y 



2 tj-i 



(Bll) 



into equation lBll The resulting differential equation takes the form 
2 



y = 



(1 + 77)70 lnA 
With the definition 

3 + 77 



ln(A j/ p ) i+„ 
V 2 



(B12) 
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which is positive for all realistic cases, we find the same implicit 
solution for (3 7^ and explicit solution for /3 — as in the Kepler 
case but with a different rat (eqs.l25landl26b. 

This paper has been typeset from a TpX/ ETpX file prepared by the 
author. 
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(B6) 



with Tdf from eq. 1281 and th e exponential-integra l function Ei(a;) 
(see 2.2724.2 and 8.211.2 of Gradshteyn, Ryzhi3 dl980h >. For the 
special case of 77 = 3/2 we can use A as variable and find with 



■ln(lnA) 



§_y 

dt"' K ' In A y 
the implicit solution 

lnAo_\ 



r d f In 



ln\(y)J 



P / 0, k = 



(B7) 



(B8) 



which can be inverted to the explicit solution given in eq. [26] 

For a point-like body the maximum and the minimum impact 
parameter D r and ago are both linear in y leading to a constant 
Coulomb logarithm (i.e. /3 = 0). For this case direct integration of 
eq, IB3l leads to the explicit solutions given in eq.1261 



